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ABSTRACT. We have determined the ice mass evolution of the Antarctica and Greenland ice sheets (AIS 
and GIS) and Gulf of Alaska (GOA) glaciers from a new GRACE global solution of equal-area surface 
mass concentration parcels (mascons) in equivalent height of water. The mascons were estimated 
directly from the reduction of the inter-satellite K-band range-rate (KBRR) observations, taking into 
account the full noise covariance, and formally iterating the solution. The new solution increases signal 
recovery while reducing the GRACE KBRR observation residuals. The mascons were estimated with 
10 day and 1 arc degree equal-area sampling, applying anisotropic constraints. An ensemble empirical 
mode decomposition adaptive filter was applied to the mascon time series to compute annual mass 
balances. The details and causes of the spatial and temporal variability of the land-ice regions studied 
are discussed. The estimated mass trend over the total GIS, AIS and GOA glaciers for the time period 1 
December 2003 to 1 December 2010 is -380 ±31 Gta -1 , equivalent to -1.05 ±0.09 mm a -1 sea-level 
rise. Over the same time period we estimate the mass acceleration to be -41 ± 27 Gta 2 , equivalent to a 
-0.11 ±0.08 mm a 2 rate of change in sea level. The trends and accelerations are dependent on 
significant seasonal and annual balance anomalies. 


INTRODUCTION 

The land-ice component of the cryosphere has important 
connections to the global climate system. Changes in land- 
ice cover impact the Earth's radiation budget through 
variations in albedo, and significantly impact the hydro- 
logical cycle as changes in water storage. The Earth's large 
areas of land ice, including the ice sheets, glaciers and ice 
caps, have the potential to significantly impact eustatic sea 
level, as they hold ~77% of the world's fresh water, which 
amounts to ~65 m of global sea-level equivalent (Barry and 
Gan, 2011). Changes in freshwater input from calving and 
melting Arctic land and sea ice can increase the haline 
forcing in the North Atlantic and impact global ocean 
circulation and climate patterns (Hakkinen and Rhines, 
2004). Variations in land-ice volume are closely correlated 
to global temperatures, as paleoclimate records show all the 
major Greenland ice sheet (GIS) changes have been correl- 
ated with temperature changes (Alley and others, 201 0). It is 
important that we accurately quantify the current spatial and 
temporal mass changes of the Earth's land ice to improve our 
understanding, modeling and prediction of its response and 
contribution to climate change. Since its launch in March 
2002, the Gravity Recovery and Climate Experiment 
(GRACE) mission has acquired ultra-precise inter-satellite 
K-band range and range-rate (KBRR) measurements between 
two co-orbiting satellites in a 450 km altitude polar orbit, 
~220km apart (Tapley and others, 2004). These data have 
vastly improved our knowledge of the Earth's time-variable 
gravity field. In particular the GRACE mission has been 
instrumental in quantifying recent land-ice mass evolution 
(e.g. Luthcke and others, 2006a, 2008; Velicogna, 2009; 


Chen and others, 201 1; Schrama and Wouters, 201 1; Jacob 
and others, 2012; King and others, 2012). 

There are numerous challenges in applying GRACE to the 
study of glacier and ice-sheet mass balances. GRACE senses 
all sources of mass variability at the Earth's surface, requiring 
removal of non-glacier changes through incorporation of 
independent datasets and models. Mass signals associated 
with glacial isostatic adjustments (GIA) and terrestrial water 
storage (TWS) are poorly known in many regions and often 
dominate GRACE error budgets. Owing to the GRACE 
mission's limited longitudinal sampling at monthly time 
intervals, filtering techniques are often applied to mitigate 
striping errors and to focus signal into a region of interest, 
often resulting in signal attenuation. When mass is 
constrained to a particular spatial extent, signal can either 
leak into or out of that domain and cause additional error. 
Even when the choice of filtering and models used to correct 
atmosphere and ocean mass variations is fixed, differences 
in GRACE solutions can occur, due to different methods 
used in processing the fundamental GRACE observations 
(Pritchard and others, 2010). 

Here we present a new study to quantify the ice mass 
evolution of the Antarctica and Greenland ice sheets (AIS 
and GIS) and Gulf of Alaska (GOA) glaciers from a GRACE 
global mascon solution. The mascon technique is a time/ 
space domain approach for describing time-variable gravity 
signals, in which time is discretized in contiguous, non- 
overlapping finite intervals and, within a given interval, 
space is divided into gravity contributions due to loading of 
mass concentrations of equivalent water (mascons) on the 
Earth surface. This study is unique and differs in many 
aspects from previous regional and global mascon solutions, 
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including those of Luthcke and others (2006a, 2008), 
Rowlands and others (2010) and Sabaka and others 
(2010). The mascons are estimated globally with 10 day 
temporal and 1 arc degree equal-area spatial sampling. Ice- 
sheet mascons are restricted to the grounded-ice areas 
contributing to eustatic sea-level change. We apply tem- 
poral and spatial constraints anisotropically, through seven 
distinct constraint regions representing the ice sheets (high- 
interior and low-elevation margins), GOA glaciers, land and 
ocean. The constraint equation is an exponential function of 
correlation time and distance applied to all possible pairs of 
mascons within a constraint region (Sabaka and others, 
2010). A detailed error analysis is provided, including 
consideration of leakage, and errors in modeling mass 
variations of the atmosphere and ocean during the GRACE 
data reduction (forward modeling). In addition, a controlled 
comparison is performed in which we use the same set of 
forward modeling and GRACE data reduction, but compute 
the mass time series using an averaging kernel filter applied 
to monthly spherical harmonic fields, as is common practice 
in the literature (e.g. Velicogna and Wahr, 2005). The 
mascon solution is iterated until convergence, and the 
impact of the iterated solution on the reduction of the 
GRACE KESRR residuals is quantified. Furthermore, this study 
differs from the mascon solution presented by Jacob and 
others (2012), in that we estimate the global mascons 
directly from the GRACE KBRR data, taking into account the 
full noise covariance, and iterate the solution while 
quantifying the minimization of the GRACE KBRR obser- 
vation residuals. We present the details of the global mascon 
solution, an error analysis and a discussion of results focused 
on recent Greenland and Antarctic ice sheet and Gulf of 
Alaska glaciers land-ice evolution. This paper also serves to 
document our mascon solution used in the international Ice 
Mass Balance Inter-comparison Exercise (IMBIE; Shepherd 
and others, 201 2). 


MASCON FORMULATION 

The GRACE processing centers provide time-variable gravity 
observations (Level-2 products) in the form of monthly sets 
of Stokes coefficients or spherical harmonic fields. These 
monthly fields are estimated directly from GRACE tracking 
data, typically without the use of constraints or a priori 
information. The use of these fields for scientific applica- 
tions, such as measuring the mass evolution of the GIS or 
hydrological basins, requires filtering due to north-south 
striping artifacts, which arise from the limited longitudinal 
sampling of the GRACE orbits within monthly time periods 
(Swenson and Wahr, 2006; Klees and others, 2008). Many 
filtering methods do not take into account the noise 
covariance of the estimated Stokes coefficients, and there- 
fore are not optimal filters (Klees and others, 2008). Loss of 
signal amplitude and limited spatial and temporal reso- 
lution, as well as signal leakage in and out of regions of 
interest, are particular problems when applying filtering 
methods. 

Instead of applying a filtering technique to the GRACE 
Level-2 monthly Stokes coefficients, we used a mascon 
technique, taking geolocatable anisotropic constraints to 
estimate the global mass change directly from GRACE KBRR 
data, which accounts for the full Stokes noise covariance. 

The formulation for mascon parameters is derived from 
the fact that a change in the gravitational potential caused by 


adding a small uniform layer of mass over a region at an 
epoch t can be represented as a set of (differential) potential 
coefficients which can be added to the mean field. The delta 
coefficients can be computed as (Chao and others, 1987): 
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where k\ is the loading Love number of degree /, to account 
for the Earth's elastic yielding which, in general, counteracts 
the additional surface density (Farrell, 1972); A, c j> and R are 
geographic latitude, geographic longitude and the Earth's 
mean semi-major axis, respectively; M is the Earth's mass; 
H(t) is equivalent water height in centimeters of the mass of 
the layer over a unit of surface area at the epoch t; ft is the 
solid angle surface area of the mascon region where /-/(f) is 
applied, dfl = cosAdAd<(>. 

The water height, H(t) (cm), can be equated to a surface 
mass density, <r(f) (kgrrr 2 ), by noting that water has a 
volume density of 1000 kg m -3 , which gives er(f) = 
10 x /“/(f). The estimated mascon parameter, Hy(f), for each 
mascon region, j, is a scale factor on the set of differential 
Stokes coefficients for that mascon region, giving a surface- 
mass delta[AUTHOR: Would it be better to say 'surface- 
mass change'?)] (cmw.e.). In matrix notation, an assemblage 
over j of the above equations provides a set of partial 
derivatives, L, of the differential Stokes coefficients, As, 
(complete to degree and order 60) with respect to equivalent 
water height, Ah, such that 

As = LAh. (3) 

The mascon parameters are part of nonlinear functions of 
the GRACE KBRR measurements and are therefore estimated 
using the Gauss-Newton (GN) method (Seber and Wild, 
1989) of nonlinear least-squares. This method seeks a series 
of differential corrections to the parameters via a lineariza- 
tion of the observation equations at the current state. This 
means that the changes in KBRR observations at the current 
state due to changes in mascon parameters are needed for 
the estimation. The partial derivatives of the GRACE KBRR 
observations with respect to the mascon parameters are 
computed by the chain rule, using the partial derivatives of 
the KBRR observations with respect to standard Stokes 
coefficients and the differential Stokes coefficients com- 
puted from the equation above for the mascon region of 
interest: 


d0i do, QA C j lm (t) dO, dAsi m (t) 

dHj(t) f^^ 0 dCi m dHj(t) dS lm dHj(t) ' 

where dO,/dHj(t) is the partial derivative of KBRR obser- 
vation / with respect to the j mascon parameter, Hj(t), at 
time t ; dO;/dCi m and dO//<9S/ m are the partial derivatives of 
the KBRR observations with respect to the geopotential 
Stokes coefficients (computed as part of the KBRR reduction 
and Level-1 data processing); d/S.C\ m (t)/dHj(t ) and 

d/S.S\ m (t)/dHj(i) are the partial derivatives of the delta 
Stokes coefficients with respect to the mascon parameter in 
region j at time f. 

In matrix notation, assemblages over / and j of the above 
equation provide a set of partial derivatives of the KBRR 
observations, o, with respect to the mascon parameters, h, 
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by chaining L with partial derivatives, A, of the KBRR 
observations with respect to the Stokes coefficients. 

Following Sabaka and others (2010), at step k of the CN 
estimation let the GRACE KBRR residuals, or, be defined as 
the difference between the observations, o, and a model 
prediction, a(x^), which is a function of a current state of a 
general parameter set, x*. In addition to the current mascon 
parameter state, this parameter set represents the effects of 
orbital arc parameters (initial baseline state orbit and 
accelerometer calibration parameters) and the mean field 
and forward models of other geophysical processes of mass 
redistribution (e.g. atmosphere and ocean tides). The details 
of the orbit arc parameters and forward models are discussed 
below, in the data processing section. The correction, Ahj,, 

to the current state of mascon parameters, h k , is then defined 
as the unique minimizer of the quadratic cost function 

J(Ah k ) =(r - ALAh^) T W(r - ALAh k ) 

+ £ i ^hj. + Ah^ Phh (h*, + Ahj. j , 

where the first term is the weighted misfit and the second 
term is a measure of mascon complexity, which in this case 
is a smoothing function that tries to drive the mascon state to 
zero. The data weight matrix, W, accounts for both 
measurement noise and the effects of the orbital arc 
parameters (back substituted), and P/,/, is a mascon 
regularization matrix, whose influence is controlled by the 

damping parameter, y,. The updated mascon state, h^ +1 , is 
given by 

{ Ah* = (L T A T WAL+^P M ) M 
GN iteration k<j (\ T A T Wr - (6) 

{ h fc+1 = h k + Ah k 

assuming a step size of unity. 

The cost function of Eqn (5) suggests a major philosophi- 
cal difference between this study and past studies, in that 
multiple iterations of GN are actually carried out. For 
instance, Rowlands and others (201 0) and Sabaka and others 
(2010) carry out a single step of GN, in which fq is zero in 
Eqn (5). In this study, the step k update, X|< +1 , which reflects 
h k+ i, is fed back into the KBRR processing in order to carry 
out the next iteration. In this sense, the minimum of 7(h) is 
sought instead of the local quadratic approximation to it at 
h kl i.e. y(Ahjt). It will be shown below that this leads to a 
significant improvement in the recovered mascon solutions, 
particularly in the ability to restore more signal strength. 

There are several ways to rewrite the GN iteration k that 
illuminate different properties of the estimator. For instance, 
the filtering nature of the estimator and its error sources can 
be seen by first expressing Eqn (6) as 

{ hjt+i = 

,TT , — i T T f r \ 

(L t A t WAL + nP hh ) L T A T w(r + ALh t J 

( 7 ) 

and then assuming the residuals can be expressed as a first- 
order linear perturbation, Ah t , about a true mascon state, 
= h^, with additive noise, e, such that r = ALAh^ + e. 
This gives 

h*+i = Rh k + Ke, (7) 


where K = (L T A T WAL + /tP hh y' L T A T W and R = KAL. The 
estimate, hj. +1 , differs from the true state, hj. +1 = + Ah<., 

due to two effects: (1) the degradation of resolving power in 
the filter, due to the presence of regularization, which is 
manifested as a deviation of the resolution matrix, R, from 
an identity matrix in the first term of Eqn (8), and (2) the 
presence of additive noise from the second term in Eqn (8). 
The effects of the resolution will be discussed in detail 
below. 

Another way to rewrite the GN iteration k is: 


s k = Lh^ 


GN iteration k< 


As k = (A t WA) Vr 
sjt+r = s k + As*, 


h, +1 =p w ;l t [/x(a t wa)- 1 +lp^l t ] 


Sfc+1 ■ 


( 9 ) 


This is now in the form of an anisotropic, non-symmetric 
(ANS) filter (Kusche, 2007) that acts upon an unfiltered 
Stokes state, s* + i , to produce a filtered Stokes state, $ k+ i, 

s*+i = p; 1 (n -1 + p *TVi, (io) 

followed by a least-square co-location (LSC) prediction 
(Moritz, 1980) of the mascon state, hjt+i , from the filtered 
Stokes state 


hfc+1 — Pfe 1 PssSit+1 / (11) 

where N _1 and P^. 1 are the noise covariance and signal 
auto-covariance matrices of the Stokes coefficients, respect- 
ively, and P^ 1 is the signal cross-covariance matrix between 
the mascons and Stokes coefficients. The ability to process 
the KBRR tracking data directly allows for access to the 
actual noise covariance matrix of the Stokes coefficients, 
and thus, makes this an example of an optimal ANS filter 
(Kusche, 2007). Note that the Stokes signal covariances are 
related to the mascon signal auto-covariance through simple 
propagation, such that P^ 1 = LP kk L T and Pj^ 1 = P kk L T . 
Currently, the mascon signal auto-covariance matrix reflects 
a first-difference operation between all distinct pairs of 
mascons in space and time, weighted as an exponential 
function of temporal and spatial distance between them, and 
a conservation of mass constraint that keeps the global 
mascon sum constant over time. Because of the regional 
aspect of the spatial weighting, the gravity signal is treated as 
a non-stationary, anisotropic process based upon geoloca- 
table physical properties. Details of the constraint regions, 
mascon signal auto-covariance construction and selection of 
the damping parameter are given in the following section. 


GRACE DATA PROCESSING, AND MASCON 

SOLUTION DETAILS 

Data processing and forward models 

We estimated global surface mass anomalies, or mascons, 
directly from the reduction of the GRACE Level-1 B data, 
which includes the inter-satellite K-band range-rate (KBRR) 
data, as well as orbit position, attitude and accelerometer 
data for each GRACE satellite. We employed a baseline state 
parameterization (Rowlands and others, 2002) and accel- 
erometer calibration (Luthcke and others, 2006b), which 
enables the estimation of surface mass variability from 
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Table 1 . Summary of forward modeling applied for each global 
mascon solution 


Mascon solution 

Forward modeling used in computation of KBRR 
residuals number and partial derivatives for 
indicated mascon solution 

v08 

GCM02C (150 x 150) 

Pole tide (IERS2003) 

Solid Earth tides (IERS2003) 

GOT4.7 ocean tides (Ray, 1999; Ray and Ponte, 
2003) 

OMCT ocean model from AOD1 B 
ECMWF 3 hour operational pressure grids 
(Nmax = 90) 

Post-solution geocenter correction (Swenson and 
others, 2008) 

v09 

Same as v08 

GLDAS/Noah 0.25°, 3 hour data (N max = 90) Ice 
'masked out' 

ICE-5G (VM2) (model C 21 and S 21 from GIA/GRACE 
solution)* 

vIO 

Same as v09 plus v09 global mascon solution 

vll 

Same as v09 plus vIO global mascon solution 

vl 2 

Same as v09 plus vl 1 global mascon solution 


*This GIA model is restored in the mascon solutions and more recent GIA 
models noted in the text are used in the computation of the mascon results. 


GRACE KBRR data alone. The GPS data are used solely to 
establish an accurate orbital reference and for calibrating 
accelerometers (Luthcke and others, 2006b). The GRACE 
KBRR and Level-1 B attitude and accelerometer data (in- 
cluding accelerometer calibration) are processed in daily 
arcs as outlined by Luthcke and others (2006b). Forward 
models for both static and time-varying gravity are applied in 
the KBRR data reduction daily arcs, in order to isolate land- 
ice signal. Static mean gravity is modeled using the 
complete GGM02C GRACE gravity model through degree 
and order (d/o) 150 (Tapley and others, 2004). Ocean tides 
are modeled with the GOT4.7 tide model (Ray, 1999; Ray 
and Ponte, 2003). Atmospheric mass variations are modeled 
to d/o 90 at 3 hour intervals using European Centre for 
Medium-Range Weather Forecasts (ECMWF) operational 
pressure grids. Ocean mass variations are modeled to d/o 90 
from the baroclinic Ocean Model for Circulation and Tides 
(OMCT). Terrestrial water storage (TWS) variations are 
modeled using the Global Land Data Assimilation System 
(GLDAS)/Noah Land Surface Model (Noah) 0.25°, 3 hour 
data, represented as potential coefficients to d/o 90 (Ek and 
others, 2003; Rodell and others, 2004). The TWS contri- 
bution was set to zero for the major mountain glacier areas 
of the globe using a 0.25° mask (Raup and others, 2000) and 
also set to zero for the GIS and AIS. Glacial isostatic 
adjustments are modeled based on ICE-5G (VM2) (Peltier, 
2004) and are represented as potential coefficients to d/o 90. 
Table 1 provides a summary of the modeling applied in the 
reduction of the KBRR residuals for each of the mascon 
solutions noted. 

GIA and geocenter 

We modified the mascon solution to restore the forward 
modeled global GIA, while more recent GIA models are 
then used in the computation of the mascon solution results 
presented. For both the GIS and GOA, a GIA model based 


on the ICE-5G deglaciation history and an incompressible 
two-layer approximation to the VM2 viscosity profile are 
used to correct the GRACE mascon solution (Paulson and 
others, 2007, computed and provided by A. Geruo, 
[AUTHOR: Please provide date]). For the AIS, the IJ05 R2 
regional GIA model is used (Ivins and others, in press). We 
correct our mascon solutions for the viscous component of 
post-Little Ice Age (LIA) GIA following the collapse of the 
Glacier Bay Icefield. As Luthcke and others (2008), we apply 
a regional post-LIA GIA model that is rigidly constrained by 
~100 precise GPS and relative sea-level (RSL) observations 
of uplift (Larsen and others, 2005). We further correct our 
mascon solutions to reflect surface mass variability in the 
Earth's center-of-figure frame, rather than the center-of-mass 
frame in which the solutions are performed. We apply the 
geocenter correction used in the IMBIE study (Shepherd and 
others, 2012) derived from the degree 1 Stokes coefficients 
determined from Swenson and others (2008). The geocenter 
correction accounts for % of the GIS and GOA ice trends, 
and <3% for the West AIS and AIS peninsula ice trends. The 
largest impact of the geocenter correction is for East AIS, at 
1 7% of the trend. 

Global mascon definitions 

We estimate a global set of 1 x 1 arc degree equal-area 
(~1 2 390 km 2 ) mascons at a 10 day temporal sampling. The 
choice of the spatial sampling is a compromise between 
approximating regional shapes and boundaries (e.g. land, 
ocean, land ice) and the total number of parameters that can 
be reasonably iterated for a global solution. There are a total 
of 41 168 mascons estimated every 10 days, representing a 
significantly smaller spatial sampling than the fundamental 
spatial resolution of GRACE (Luthcke and others, 2006b). 
Furthermore, the number of mascon parameters is signifi- 
cantly larger than the number of d/o 60 Stokes coefficients 
forming our basis functions (3718). Therefore, as discussed 
in detail above, we employ a regularization in the form of 
the mascon signal auto-covariance matrix (Sabaka and 
others, 2010). This matrix is constructed from constraint 
equations which require all mascon differences to be close 
to zero in a statistical sense. The constraints are applied in 
groups that pertain to geographical regions, so if two 
mascons reside within the same region, the constraint 
weight is a simple exponential function of the distance 
(characteristic length) and time (characteristic time) between 
the two mascons. If two mascons reside in different regions, 
the weight is zero and the constraint vanishes. The details of 
the constraints and the construction of the regularization 
matrix, P/,/,, are given by Sabaka and others (2010, their 
section 2.2). The important implementation parameters are 
the characteristic length and time used in the exponential 
taper weighting of the mascon first-difference constraint 
equations, and the damping parameter, which controls the 
influence of the regularization relative to the data misfit. 

Application of regional constraints 

The constraints are then applied anisotropically using the 
following seven constraint regions: (1) GIS margins defined 
as <2000 m elevation, (2) GIS >2000 m elevation, (3) AIS 
<2000 m elevation, (4) AIS >2000 m elevation, (5) GOA 
glacier area, (6) land (including all other ice regions) and 
(7) ocean and large bodies of water (including ice shelves). 
In defining the mascons for land and large bodies of water, a 
minimum enclosed region area of 80 000 km 2 was used. 
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Fig. 1 . Land-ice mascon configurations and drainage systems, (a) The CIS. The numeric labels describe the primary drainage systems and 
subregions with the first and second digits respectively, and follow definitions based on those of Zwally and Giovinetto (201 1). The yellow 
border delineates the 2000 m elevation cutoff used in constraining the mascon solution, (b) The AIS. The drainage systems are sequentially 
labeled from 1 to 36. The 2000 m elevation line defining the constraint regions is again delineated by a yellow border. The three Antarctic 
regions of EAIS (drainage systems 3-15, 26-36), WAIS (drainage systems 1, 2, 16-22) and AIS peninsula (drainage systems 23-25) are 
outlined in red. (Note that the yellow elevation line is concealed by the red line along a portion of the boundary between EAIS and WAIS.) 
(c) The GOA mascons. This region is not further delineated by drainage systems or constraint regions. 


Thus, if a land area enclosed within the ocean region is 
<80 000 km 2 , then those land mascons were labeled as 
ocean (similarly for water regions enclosed within land). 
Figure 1 shows the GIS, AIS and GOA mascons as well as 
the constraint region boundaries. Also shown in this figure 
are ice-sheet drainage system divides based on those of 
Zwally and Giovinetto (2011). These divides are used for 
computing mass change as the sum of the mascon time 
series within a drainage system, but are not used in the 
mascon estimation. 

We use a characteristic length of 1 00 km and a character- 
istic time of 10 days in computing the exponential 
constraints for the mascon signal auto-covariance matrix. 


These values represent approximately one parameter sample 
in time and space. We use a damping parameter {fi in 
Eqn (5)) of 30{,000~\ which provides a similar solution 
contribution from the signal and data covariance, as 
determined from an analysis of the GIS, AIS and GOA 
mascon matrix diagonals. We investigated the impact of the 
damping parameter by performing solutions using values 
that span an order of magnitude variation, centered on our 
selected damping parameter value. We find that the impact 
on the estimated mass balance and annual amplitude is 
<8% on average over the GIS, AIS and GOA mascons, while 
the 1 G error estimates of the 10 day mascon solutions can 
change by as much as 31% on average over those regions. 
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Fig. 2. Daily KBRR residual rms solution differences. The daily 
(dots) and 15 day Gaussian smoothed (curves) of the KBRR residual 
rms differences between the v08 solution and the v09 (blue), vIO 
(purple), vll (green) and v12 (orange) solutions. Negative values 
represent a reduction in the residuals and an improvement in the 
forward modeling, including an iteration of the mascon solution. 

Therefore, varying the signal covariance damping parameter 
has a larger impact on the solution noise than the signal. 

Resolution and leakage 

We further analyze the impact of our signal covariance 
(including our choice of the damping parameter and 
characteristic length and time) by performing signal leakage 
analysis within and across constraint regions using the 
resolution matrix identified in Eqn (8). First, we look at the 
influence of the signal covariance on individual mascons. 
Several GIS and AIS mascons were chosen as single-source 
mascons. For each single-source test we fill the mascon 
parameter vector with a time series of mass change for that 
mascon only (we use the time series from our final solution) 
and zero all other mascons in space and time. The resolution 
matrix then multiplies the vector of mascons in time and 
space, where only one mascon contains a time series of 
mass change. Therefore, before application of the resolution 
matrix, only one mascon out of all the global mascons 
contains a nonzero time series of mass change. However, 
after application of the resolution matrix several neighboring 
mascons will contain nonzero time series, as the signal 
covariance spreads out the signal from the source mascon. 
The effects of the signal covariance are quantified by 
differencing the full mascon vector before and after the 
application of the resolution matrix. The results of several 
tests show that the spatial signal from a single mascon is 
~100% accounted for by adding all the mascons within a 
600 km radius of the original source mascon. Therefore, a 
single mascon has negligible influence at distances 
>600 km. This is equivalent to a Gaussian spatial smoothing 
filter with a 300 km radius. However, the constraint region 
boundaries modify the spatial pattern of signal spreading as 
the signal covariance, by design, limits signal leakage across 
these boundaries. The spatial leakage pattern of a single 
mascon's signal is a function of its location with respect to 
the constraint region boundaries. An error analysis of signal 
leakage across constraint boundaries is performed for each 


region considered. The regional leakage analysis and its 
results are discussed further below. In summary, the spatial 
resolution of our solutions is equivalent to a 300 km 
Gaussian spatial smoothing filter within constraint regions, 
but with significantly reduced leakage across constraint 
region boundaries. 

MASCON SOLUTION ITERATION AND KBRR 
RESIDUAL ANALYSIS 

As described above, we apply a mascon technique to 
estimate global mass change directly from the GRACE KBRR 
observations themselves, taking into account the full Stokes 
noise covariance and a signal covariance. Therefore, we can 
use the reduction in the GRACE KBRR residuals as a 
quantitative objective measure of performance. The KBRR 
residuals are the difference between the GRACE observed 
KBRR data and KBRR data computed from our forward 
models and mascon solution using the NASA Goddard 
Space Flight Center GEODYN precision orbit determination 
and geodetic parameter estimation software. Unique to this 
study, we present the improvement in the KBRR residuals 
from fully iterating the global mascon solution (as in 
Eqn (6)). Because there are > 1 1 x 1 0 6 mascon parameters, 
we obtain the solution using iterative methods. Following 
Sabaka and others (2010), we use the pre-conditioned 
conjugate gradient method. Figure 2 presents the KBRR 
residual performance for each solution defined in Table 1. 
This figure shows the difference in the daily root-mean 
square (rms) of the KBRR residuals using the v08 forward 
model and the KBRR residuals computed from subsequent 
versions of forward models (which include iterated mascon 
solutions for v10-v12). Negative differences represent 
improvements due to the new forward model of mass 
change. The v09-v08 daily rms difference in Figure 2 
demonstrates significant reduction in the KBRR residuals 
(model improvement) from forward modeling hydrology and 
GIA. The residual reduction is seasonal, with the largest 
reduction in daily residual rms in Northern Hemisphere 
spring and an increase in residuals during a portion of 
Northern Hemisphere summer. The residuals show that the 
GLDAS hydrology model improves the global surface mass 
modeling during most of the year, except the Northern 
Hemisphere summer. 

The global mascon solutions represent a reduction in the 
daily KBRR residual rms, as shown in Figure 2. The first 
iteration mascon solution, vIO forward model, represents a 
factor of 2.4 reduction in residuals over the v09 forward 
model (v08 plus GLDAS hydrology and GIA), and further 
iterations represent 5% (vll forward model) and 2 % (v12 
forward model) reduction over the preceding iteration. The 
global mascon solutions show reduction in seasonal and 
long-term trend signals in the daily KBRR residuals, due to 
improved modeling of the global mass change, dominated 
by hydrology and land ice. The spatial distribution of the 
residual reduction is presented in Figure 3, as the trend and 
annual amplitude of the GRACE average time-differenced 
KBRR residuals (KBRR-dot). The significant trends in the 
KBRR-dot residuals correspond to the GOA, AIS and GIS 
land-ice regions, and in particular, the largest residual trends 
are found in the West Antarctic ice sheet (WAIS) Amundsen 
and Bellingshausen coasts, and the southeast and northwest 
coasts of the GIS (Fig. 3a). Figure 3a demonstrates the 
significant reduction in these KBRR-dot residual trends using 
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Fig. 3. Global maps of the (a) trend and (b) annual amplitude of the time-differenced KBRR residuals (KBRR-dot) for the forward models used 
in the v08 (row 1), v09 (row 2) and v12 (row 3) mascon solutions (Table 1). The differences between v09 and v08 show the impact of the 
additional forward models included in the v09 solution, while the differences between v12 and v09 show the large reduction in residuals 
that is achieved with the iterated mascon solution. 


the v12 forward model, due to the inclusion of the vl 1 
iterated global mascon solution. 

Iterating the solution has an important impact on the 
regional terrestrial ice mass signals obtained from the global 
mascon solution (Fig. 4 and Table 2). Iteration results in 
increases in ice mass loss trend, annual and acceleration 
signal computed from the global mascon solution. It 
accounts for a —16.5, —21.8 and —1.1 Gta -1 correction to 
the GRACE observed GIS, AIS and GOA mascon solution 
trends; and 1 0.7, 34.2 and 4.4 Gt corrections to the GIS, AIS 
and GOA annual amplitude (Table 2). Iteration of the global 
mascon solution is shown to reduce KBRR residuals and 
increase signal recovery. Therefore, we regard the increased 
signal recovery through iteration as an important part of the 
global mascon solution, without which the regional land-ice 
mass signals estimated from the mascon solution would be 
in error by the amounts noted in Table 2 (v12-v09). 

SIGNAL AND ERROR ANALYSIS 

For each of the 41 168 mascons in our global solution, we 
obtain a time series of mass anomalies with 10 day 
sampling. To compute the time series of mass anomalies 
for any region of interest we sum the individual time series at 
each 10 day sample for all mascons contained within the 


region. The 10 day mascon solution Icr error is estimated 
from the rms of the misfit of the mascon solution time series 
and a 1 0 day width Gaussian filter applied to the time series, 
in a similar way to Luthcke and others (2008). Also, as in 
Luthcke and others (2008), the mascon time series trend and 
annual amplitude are determined from a least-squares 
simultaneous estimate of the trend, annual, semi-annual 
and 1 61 day periodic. The 1 61 day periodic compensates for 
the GRACE alias of an S 2 tide error (Ray and Luthcke, 2006). 
For the trend, we replace the formal error estimate (variance 
of the least-squares parameter) with one that considers a 
first-order autoregression structure for the mascon time 
series (Lee and Lund, 2004). The formal estimates of the 
other parameters (e.g. annual amplitude) are then scaled by 
the ratio of the trend error with autoregressive structure 
considered, and the simple formal trend error. The total error 
is then computed as the root-sum-square (rss) of the 
following error estimate sources: statistical; leakage in and 
out; atmosphere/ocean forward modeling and GIA. These 
errors are discussed below. 

Leakage errors 

Leakage of signals in and out of regions of interest occur clue 
to the fundamental spatial resolution limitations of the 
GRACE data itself and the properties of our mascon signal 
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Fig. 4. The effect of iteration on the mass time series for (a) CIS, (b) AIS and (c) GOA. The third and final iteration solution is shown (v12: 
gray), along with the differences to the original (v09: blue), the first iteration (vl 0: purple) and the second iteration (vl 1 : green) solutions. 


autocovariance. Errors arising from leakage of signals into 
('leakage-in') a region of interest (e.g. GIS or AIS) from 
unmodeled surrounding mass variation sources, such as 
ocean, tides and hydrology, are estimated, as well as land- 
ice mass signal leakage out ('leakage-out') of our regions of 
interest. It should be noted that the hydrology signal leakage 
error has been significantly reduced through the forward 
modeling of a TWS model, as previously discussed, and as 
in Luthcke and others (2008) and Sabaka and others (2010). 
We use a representative global mass anomaly model in 
space and time to estimate the leakage-in and leakage-out 
errors. For this model, we use the vector (Eqn (8)) from 
our final perferred v12 global mascon solution. To estimate 
leakage-in for a region of interest, we zero the region 
mascon parameters in hj. from Eqn (8), but leave the 
remaining global mascon parameter values to produce 
h/y?(o). The estimate of leakage-in error is then the region's 
mascon time series after applying the resolution operator to 
our model, h^o). We then perform this procedure for each 
region of interest. The leakage-out errors are computed using 
a similar procedure, but in this case all mascon parameters 
outside the region of interest are set to zero. The estimate of 
leakage-out error is then the difference between the region's 
model mascon parameters and those after applying the 
resolution operator. A summary of the leakage errors for the 
major land-ice regions of interest is given in Figure 5a and 
Table 3. While the leakage analysis does consider sources 
(1 arc degree mascons) significantly smaller than the 
fundamental spatial resolution of GRACE (approximately 
equivalent to 300 km radius Gaussian smoother), the analy- 
sis did not consider point sources and localized signal 
smaller in spatial extent than an individual mascon. An 
assessment of the impact of sub-mascon sources on our 
leakage estimates will need to be evaluated in future 
analysis. 

Atmosphere/ocean forward modeling errors 

Following Han and others (2004), we estimate the errors in 
forward modeling atmosphere and ocean mass as the 
difference between two sets of models: (1) ECMWF and 
OMCT vs (2) atmospheric mass variations derived from 
NASA/GSFC Modern Era Retrospective-analysis for Research 
and Applications (MERRA) model and a simple ocean 
inverted barometer (IB). A summary of the impact of these 
models on our land-ice solutions is provided in Figure 5b 


and Table 4. The two models are not completely inde- 
pendent, and, therefore, we do not scale the differences by 
1/y/2 (Han and others, 2004). Ocean mass variation and 
ocean tides forward modeling error contributions to our 
land-ice mass solutions are included through the leakage-in 
analysis discussed above. The mass anomaly model used in 
the leakage analysis includes the ocean mass and ocean 
tides modeling error observed by GRACE, as the mascon 
solution reflects the mass anomalies from the forward 
models. 

GIA errors 

The land-ice mass trends estimated from our global mascon 
solutions contain errors due to the uncertainties in the GIA 
model applied. We estimate the GIA error as the difference 
between the minimum and maximum GIA correction over a 
range of plausible GIA models divided by 2 for each region 
of interest. For the GIS and GOA the models considered 
include: the nominal ICE-5G Paulson incompressible 2-layer 


Table 2. Effect of solution iteration on mascon regional mass change 
statistics over the time period 1 December 2003-1 December 
2010. The reported statistics describe the differences of the final 
and third iteration solution (vl 2) to the original global solution 
(v09), the first iteration solution (vIO), and the second iteration 
solution (vl 1) 


Solution iteration 
by region 

Annual 

amplitude 

Gt 

Acceleration 

Gta -2 

Trend 

Gta- 1 

CIS 

vl 2-v09 

10.7 

-0.7 

-16.5 

v12-v1 0 

6.3 

-1.0 

-7.8 

vl 2-vl 1 

2.8 

-0.5 

-3.2 

AIS 

vl 2-v09 

34.2 

-11.3 

-21.8 

vl 2-vl 0 

18.3 

-5.8 

-10.3 

vl 2-vl 1 

7.6 

-2.5 

-4.2 

GOA 

vl 2-v09 

4.4 

2.9 

-1.1 

vl 2-vl 0 

1.9 

1.6 

-0.9 

vl 2-vl 1 

0.8 

0.6 

-0.5 
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Fig. 5. (a) The vl2 mascon solution mass time series for GIS, AIS and GOA with leakage-in (black circles) and leakage-out (red diamonds) 
estimated errors, (b) The v12 mascon solution time series for GIS, AIS and GOA with atmosphere and ocean model signal (ECMWF/OMCT: 
black circles) and estimated model errors (ECMWF/OMCT - MERRA/IB: red diamonds), (c) The GIS v12 mascon solution time series and the 
difference to the mass time series derived from the averaging kernel filter technique (red diamonds). 


model discussed previously; an ICE-5G compressible model 
(Geruo and others, in press) and models based on 'ANU' ice 
deglaciation history (Fleming and Lambeck, 2004) applying 
a suite of Earth models (27) over a range of lithospheric 
thicknesses, and upper mantle and lower mantle viscosities 
(Barletta and others, 2012, computed and provided by 
V.R. Barletta). For the AIS the models considered include: 
the nominal IJ05 R2 used in this study and discussed 
previously; models based on IJ05 R2 with varying Earth 
models and models based on W12a with the 'best'-fit Earth 
model and upper- and lower-bound Earth models that give a 
good fit to the relative sea-level data (Whitehouse and 
others, 2012, computed and provided by P. Whitehouse). 
This is a similar approach to estimate the GIA uncertainty as 
that used in the IMBIE study (Shepherd and others, 2012). 
The post-LIA GIA signal uncertainty is ±30%, estimated 
from the model comparisons to GPS vertical motions, raised 
shoreline records of RSL and tide-gauge rates of RSL (Larsen 
and others, 2005). 

Solution technique error analysis 

As another step in our solution evaluation, we compare our 
mascon solution with that of a solution made using the more 


common approach of applying a filter to the monthly 
spherical harmonic fields. In this test we apply a calibrated 
regional GIS averaging kernel to our monthly spherical 
harmonic fields estimated from the same Level-1 B proces- 
sing as our mascon solution (Luthcke and others, 2006b). 
The averaging kernel is constructed, calibrated and applied 
in a similar fashion to Velicogna and Wahr (2005). The 
important aspect of this comparison experiment is that both 
the mascons and the monthly spherical harmonic fields have 
been estimated from the same Level-1 B processing and 
KBRR reduction, and, therefore, the results isolate the 
differences arising from the technique used to compute the 
total GIS mass time series. Figure 5c shows (in red) the 
difference between the time series computed from the GIS 
mascons and that computed from the averaging kernel; the 
mascon solution time series (black) is also provided in this 
figure for comparison. The differences in the annual ampli- 
tude, acceleration and trend between the two time series 
are, respectively, 25 Gt, -1 Gta~ 2 and -lIGta -1 . These 
potential uncertainties are due solely to the differences in 
the techniques used for the mascon time series recovery, and 
are on the order of the estimated uncertainties for our 


Table 3. Effect of leakage-in, leakage-out and the rss of the leakage-in and -out, on the annual amplitude, acceleration and trend of mass 
changes for each main region and subregion over the time period 1 December 2003-1 December 2010 


Region 

Annual 

amplitude 

Gt 

Leakage- in 
Acceleration 

Gta -2 

Trend 

Gta -1 

Annual 

amplitude 

Gt 

Leakage-out 

Acceleration 

Gta- 2 

Trend 

Gta- 1 

Annual 

amplitude 

Gt 

rss leakage 
Acceleration 

Gta- 2 

Trend 

Gta- 1 

GIS total 

8.7 

-2.9 

2.5 

17.6 

-0.2 

-7.3 

19.7 

2.9 

7.7 

GIS <2000 m 

23.8 

-2.14 

-3.6 

35.8 

1.9 

-17.9 

43.0 

2.8 

18.2 

GIS >2000 m 

15.3 

3.2 

-7.1 

12.2 

1.9 

-2.6 

19.6 

3.7 

7.6 

AIS total 

12.4 

0.5 

-10.9 

16.8 

0.1 

-4.4 

20.8 

0.6 

11.8 

AIS <2000 m 

15.3 

2.9 

-3.6 

16.3 

1.4 

-2.1 

22.4 

3.2 

4.2 

AIS >2000 m 

2.8 

1.0 

-3.3 

5.9 

2.0 

1.7 

6.5 

2.2 

3.7 

AIS east 

6.7 

0.7 

-11.1 

8.6 

1.6 

1.3 

10.9 

1.8 

11.2 

AIS west 

5.1 

-1.7 

-1.7 

5.4 

-3.9 

-8.6 

7.4 

4.3 

8.8 

AIS peninsula 

3.0 

-2.7 

-10.0 

5.3 

-1.8 

-9.1 

6.1 

3.3 

13.5 

GOA total 

19.0 

1.4 

6.9 

34.3 

0.8 

4.2 

39.2 

1.6 

8.1 
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Table 4. Annual amplitude, acceleration and trend of mass changes due to the atmospheric and ocean signal and errors for each main region 
and subregion over the time period 1 December 2003-1 December 2010 


Atmospheric and ocean signal* Atmospheric and ocean errod 


Region 

Annual amplitude 
Gt 

Acceleration 

Gta -2 

Trend 

Gta- 1 

Annual amplitude 
Gt 

Acceleration 

Gta“ 2 

Trend 

Gta” 1 

GIS total 

199.1 

7.1 

3.5 

9.5 

-2.5 

-4.5 

GIS <2000 m 

90.6 

2.8 

-0.3 

8.9 

-2.0 

-4.1 

GIS >2000 m 

109.9 

4.3 

3.7 

1.7 

-0.5 

-0.4 

AIS total 

697.7 

-58.5 

-5.8 

17.6 

-21.5 

8.0 

AIS <2000 m 

262.9 

-40.4 

-16.9 

7.7 

-13.5 

-6.5 

AIS >2000 m 

437.0 

-18.1 

11.1 

10.3 

-7.9 

14.5 

AIS east 

561.2 

-37.4 

2.2 

15.1 

-17.6 

8.3 

AIS west 

127.6 

-17.7 

-5.9 

8.3 

-2.2 

-0.1 

AIS peninsula 

9.9 

-3.4 

-2.1 

4.2 

-1.7 

-0.2 

GOA total 

42.8 

-0.7 

-2.1 

1.7 

-0.6 

-1.8 


*Sum of ECMWF atmospheric model and OMCT ocean model. 

Difference between the combined ECMWF/OMCT models and the combined MERRA/IB models. 


mascon solution, as summarized in Table 5 and discussed in 
the next section. 

Computation of seasonal balances 

We use a linear regression analysis to extract the multi-year 
average trend, acceleration and annual amplitude from our 
10 day mascon time series. One of the strengths of GRACE 
solutions, and of the mascon solution presented here in 
particular, is the high temporal resolution at which the 
evolution of land-ice mass balance can be observed and 
quantified. Therefore, we seek to compute annual mass 
balances from our GRACE mascon time series, in order to 
better quantify the more complex spatial and temporal 
variability of land ice. Luthcke and others (2008) simply 
used the extrema of the full mascon time series as the start 
and end of the balance years from which to compute the 
annual mass balances. Because the mass time series for the 
sum of GOA glacier mascons has a very clean annual signal 


with little additional signal or noise, this is an acceptable 
approach. However, for the more complex mascon time 
series from the GIS, WAIS, East Antarctic ice sheet (EAIS) and 
Antarctic ice sheet peninsula (AISP), an adaptive filter 
approach is necessary. Therefore, we employ an ensemble 
empirical mode decomposition (EEMD) adaptive filter to first 
compute the timing of the balance seasons (accumulation 
and loss), and then to compute the annual mass balances. 

The EEMD is a noise-assisted data analysis improvement 
of the EMD method, which provides an adaptive decom- 
position of the time series into components known as 
intrinsic mode functions (IMFs) (Wu and Huang, 2009). The 
IMFs are monochromatic at a given point in time, but vary in 
time with both amplitude and frequency. The IMF adaptive 
decomposition separates out the frequency content of the 
time series, which represent different physical modes. An 
ensemble (many runs) of the EMD is performed, applying a 
new generation of noise to the original signal on each run 


Table 5. Table 5. Summary of the v!2 global mascon solution for the land-ice regions studied over the time period 1 December 2003-1 
December 2010. The reported parameters are computed after applying the GIA and LIA model correction to the mascon time series. The 
total uncertainties include contributions from statistical errors, mismodeling of the atmosphere, leakage errors, GIA model errors and, for 
Gulf of Alaska, the LIA model errors 


Region 

Area 
1 0 6 km 2 

Annual amplitude 
Gt 

GRACE mascon solution 
Acceleration 
Gta- 2 

Trend 

Gta -1 

GIA model* 
Trend 
Gta- 1 

GIS total 

2.46 

131 .9±28.6 

-10.1 ±6.3 

-230.3 ±12.1 

4.7 ±6.7 

GIS <2000 m 

1.49 

109.0±48.6 

-1 0.6 ± 6.7 

-223.7 ±19.8 

6.9±3.8 

GIS >2000 m 

0.98 

23.1 ±22.3 

0.5 ±4.7 

-6.6 ±8.6 

-2.2 ±3.0 

AIS total 

13.96 

31 6.4 ± 57.7 

-32.6±25.6 

-80.8 ±26.2 

47.0 ±17.9 

AIS <2000 m 

6.93 

1 72.6 ± 50.7 

-32.6 ±18.6 

-1 1 3.3 ±21 .0 

42.0 ± 1 5.9 

AIS >2000 m 

7.02 

143.7 ±25.1 

0.1 ±10.2 

32.5 ±21.6 

5.0 ± 1 4.6 

AIS east 

10.57 

215.4±48.8 

22.3 ±21 .6 

62.8±28.1 

1 6.8 ±21.6 

AIS west 

2.76 

83.2 ±17.6 

-45.5 ± 6.1 

-1 06.0 ±15.9 

26.2 ±12.8 

AIS peninsula 

0.63 

1 8.3 ±14.6 

-9.3 ±5.1 

-37.7 ±14.0 

4.1 ±1.6 

GOA totaft 

0.76 

1 65.9 ±46.2 

1. 7*6.9 

-68.8 ±11.0 

2.1 ±2.0 


*GIS and GOA are corrected with the ICE5G_Paulson GIA model. AIS is corrected with the IJ05_R2 GIA model. 
iThe trend correction for the GOA region from the Larsen LIA model is 10.2 ±3.1 Gta~h 
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Fig. 6. The v12 10 day mascon solution mass time series (circles) 
and 10 day Gaussian smoothed solution (curves) for GIS total 
(blue), GIS <2000 m elevation (purple) and GIS >2000 m elevation 
(green). The red lines show the best-fit trend plus acceleration. 

within the ensemble. The added finite white noise provides a 
means to sort out signals with similar scale into their proper 
IMFs and thus mitigate mode mixing. The mean of the 
ensemble is treated as the final true result, and, therefore, 
the effect of noise is minimized. Unlike typical signal 
analysis techniques, such as the linear regression analysis 
applied above, the EEMD is an adaptive filter and does not 
rely on linear and stationarity assumptions. 

We apply the EEMD to adaptively find the extrema of the 
non-stationary annual signal within the mascon time series of 
interest. The minima and maxima of the annual IMF represent 
the timing of the mass-balance seasons (beginning of the 
mass accumulation and mass loss seasons). This is the most 
important aspect of our EEMD application. The adaptive filter 
is used to isolate the annual signal, and then to identify the 
timing of the balance seasons and length of each balance 
year as the annual signal is non-stationary and evolving. The 
original signal is then reconstructed by summing the IMFs, 
and eliminating noise by discarding the first few IMFs before 
the IMF preceding the annual[AUTHOR: Please clarify 
meaning of 'preceding the annual'] (Fig. 8). We then 
compute the mass value of the reconstructed signal at the 
beginning of the mass-balance seasons (annual IMF minima). 
We use 100 runs within each EEMD ensemble, where the 
added white noise is recomputed on each run within the 
ensemble, with a noise ratio equal to twice the ratio of the 
standard deviation of the mascon solution noise (difference 
between 10 day solution time series values and a 30 day 
Gaussian filter), and the standard deviation of the signal (with 
trend removed). The EEMD is then run 50 times, where the 
overall signal trend is removed before each run and is added 
back after the reconstruction of the signal from the selected 
IMFs. The mean (computed over the 50 EEMD runs) timing 
value of the annual IMF minima and corresponding 
reconstructed signal mass values are used to compute the 
annual balance (Fig. 8). The annual balance is computed by 
differencing the reconstructed mass values at the times of two 
consecutive mean annual IMF minima. The mass difference, 
or annual balance, is labeled with the year that is mostly 
encompassed between the two minima. Therefore, for the 
GIS and GOA land ice, the balance year starts in the late 


summer to late fall [AUTHOR: Should this be 'early fall'?] of 

the previous year, while for the AIS the balance year starts in 
the Northern Hemisphere spring. The standard deviation of 
these values over the 50 runs is used as an estimate of the 
noise contribution to the uncertainty. To estimate a potential 
systematic error in our balance year timing and values, the 
EEMD analysis results are differenced from an iterative 
piecewise trend and annual linear regression analysis of the 
mascon time series, where the annual pieces are constrained 
to form a continuous signal. The error estimate of balance 
year extrema timing and mass value is then computed as the 
root-sum-square (rss) of the noise and systematic error 
estimates. The error estimate of the annual balance is then 
the rss of the error estimates for the two balance year 
extrema. 


LAND-ICE MASS EVOLUTION RESULTS 

In this section we provide a summary of our iterated global 
mascon solution, focusing on the results for the five most 
significant land-ice regions: GIS, WAIS, EAIS, AISP and 
GOA. Table 5 provides an overall summary of the regional 
mass time series computed from the mascon solution. While 
the time series cover a longer time period, the model 
parameter fit is applied to an integer number of years, 1 
December 2003 to 1 December 2010. We use the 
terminology defined by Cogley and others (201 1) to explain 
glacier and ice-sheet variations in terms of variations in their 
climatic balances, S c | im , and calving fluxes, D, ignoring the 
contributions from basal balances. 

GIS results 

The time series for the sum of the GIS mascons (GIS total), 
and the mascons above and below 2000 m are shown in 
Figure 6, and a summary of these time series is provided in 
Table 5. The spatial distribution of the GIS trend and 
acceleration is shown in Figure 7. The GIS exhibits the 
largest rates of mass loss relative to both Antarctica and 
Alaska, with a trend of —230 ± 12 Gta -1 , an acceleration of 
mass loss of — 10±6Gta~ 2 , and an average annual ampli- 
tude of 132±29Gt. Both the magnitude of the GIA 
correction and our estimated uncertainty determined from 
the model differences are quite small in comparison with the 
total corrected trend (Table 5). Our estimates agree to within 
stated error bars with previous estimates spanning at least 
the period 2003 to October 2009: -219±38Gta _1 for 
April 2002-November 2009 (Chen and others, 2011); 
—238 ± 29 Gta -1 for October 2003-October 2009 (Sasgen 
and others, 2012); — 222 ±9Gta~ 1 for January 2003- 
December 2010 (Jacob and others, 2012). Our estimated 
acceleration of — 10±6Gta~ 2 is in agreement with 
— 15±7Gta~ 2 for August 2001-August 2010 (Sasgen and 
others, 2012) and — 1 7 ± 8 Gta~ 2 for April 2002-June 201 0 
(Rignot and others, 2011). The GIS mass evolution is 
dominated by the signal from the coastal margins of the 
ice sheet, as defined in our solution as mascons <2000 m. 
Large coastal mass loss rates result from high rates of D 
associated with the presence of fast-flowing marine-termi- 
nating glaciers along the northwest, southwest and southeast 
regions of the GIS (Moon and others, 2012), and due to the 
higher surface melt rates that occur at low elevations on the 
ice sheet. Above 2000 m there is little statistically significant 
signal. 
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Fig. 7. The GIS v!2 mascon solution mass changes computed over 
the time period 1 December 2003-1 December 2010. (a) Trend 
corrected with ICE5G_Paulson GIA model and (b) acceleration. 

The GIS annual balance results from the application of the 
EEMD filter are shown in Figures 8 and 9. The GIS annual 
balance anomalies (difference from the mean) show large 
interannual variability (Fig. 9). The GIS mean annual balance 
over the balance years 2004-10 is — 248 ± 12 Gta~' . It is 
important to note that our mean annual balance is slightly 
different from the trend reported above, primarily due to the 
differing time periods, where integer balance years are 
needed for the annual balance computation. The balance 
years 2005-09 have anomaly magnitudes <20% of the mean 
annual balance and alternate between positive and negative 
values. In contrast, the GIS balance years 2004 and 2010 
have large positive and negative anomalies, ~50% of the 
mean annual balance, giving rise to the statistically signifi- 
cant acceleration determined from the linear regression 
analysis (Table 5). Therefore our reported GIS mass loss 
acceleration is, in part, an artifact of the GRACE sampling 
interval, rather than an indication of a persistent multi-year 
signal, and longer time series will be needed before making 
inferences about future evolution of the ice sheet. 

In the northwest GIS (drainage systems 81 and 82) we 
observe an increase in mass losses after 2005 (Figs 10 and 
1 1), attributed to a reduction in snow accumulation (Sasgen 


and others, 2012) that combines with increased D after 

2007 (Moon and others, 2012) to further enhance mass 
losses. The timing and extent of the northward progression of 
mass loss is confirmed by GPS observations, showing 
enhanced rates of uplift late in 2005 at Thule Station in 
northwest Greenland (Khan and others, 2010). In the 
southwest GIS (drainage systems 61, 62 and 72) we observe 
persistent mass losses, attributed to high D that commenced 
in the late 1990s, that are dominated by the large calving 
losses from Jakobshavn Isbrae, the largest marine-terminating 
outlet on the GIS Qoughin and others, 2012). 

Large interannual mass variability is observed in drainage 
system 32, the region of the Geikie Peninsula in central East 
Greenland that is largely composed of ice not receiving flow 
contributions from the GIS (Fig. 10) Oiskoot and others, 
2012). The mass time series for this region shows an increase 
in mass loss beginning in 2005, followed by mass gains from 
2007-09 (Fig. 11). The glaciers in this area are smaller and 
steeper than ice-sheet outlet glaciers and exist in a region 
with high precipitation rates (Ettema and others, 2009), so 
they likely respond rapidly to climate variations in a fashion 
similar to other maritime glacier systems. Approximately 
90% of the glaciers in this region drain through tidewater 
glacier outlets (Jiskoot and others, 2012), and those draining 
to the south into the Blosseville Kyst may be particularly 
sensitive to changes in climate, due to their southerly aspect 
and concentration of ice at lower elevations. 

Several studies have reported synchronous retreat of 
marine-terminating glaciers in southeast Greenland during 
2000-06 (Howat and others, 2008), followed by advance 
and thickening of the Helheim and Kangerdlugssuaq Gla- 
ciers (Howat and others, 2007; Joughin and others, 2008). 
The reduction in mass loss associated with these advances 
has been assumed to be synchronous across the southeast 
coast (Murray and others, 2010), although evidence for this 
was based, in part, on terminus advance rates that did not 
measure glacier mass or volume changes directly (Howat 
and others, 2008). Our results show a deceleration in mass 
loss during the entire 2003-10 period (Fig. 7b), but analysis 
of the individual mass-balance seasons shows that much of 
the deceleration occurred due to mass gains in central east 
GIS in 2008, and a synchronous reduction in mass losses in 
2009 (Fig. 10). These findings are consistent with those of 
Chen and others (2011), who used GRACE observations to 
infer a significant reduction in the mass loss of southeast GIS 
glaciers during September 2007-November 2009, however 
our results pinpoint the timing of this reduction in mass loss 
to the 2009 balance year. The 2009 slowdown of mass loss 
in southeast GIS likely occurred due to a complex array of 
factors, including a deceleration and advance of several 
southeast GIS glaciers (Murray and others, 2010), as well as 

2008 and 2009 positive accumulation anomalies (Sasgen 
and others, 2012). We note that observed deceleration and 
advance of southeast GIS glaciers occurred one to two years 
prior to our observed 2009 reduction in mass losses. This 
reflects the complex and asynchronous nature of ice-sheet 
dynamics in southeast GIS (Moon and others, 2012). 

Our observations show evidence of the response of the 
GIS to variations in summer air temperatures and the length 
of the summer melt season. We observe high rates of mass 
loss across western and northwestern GIS during 2008, 
associated with extreme summer snowmelt that summer, 
detected by passive microwave imaging sensors (Tedesco 
and others, 2008). Coastal mass balances during 2009 were 
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Fig. 8. The EEMD v12 mascon solution time series and error analysis for (a) CIS, (b) GOA, (c) EAIS, (d) WAIS and (e) AIS peninsula. Each plot 
shows the v12 mascon solution time series (black circles), the mean of 50 EEMD runs, where each case sums the IMFs from the one 
preceding the 'annual' IMF to the final IMF (blue curve), the mean values from the 50 EEMD runs at the times of the 'annual' IMF extrema 
(red triangles) and the extrema error ellipses (green ellipses). 


less negative than all other years in the GRACE record, and 
this was followed in 2010 with the most negative mass 
balances <2000 m (Fig. 10). These patterns agree well with 
observed and simulated melt extent records showing low/ 
high melt extents in the 2009/10 summers (Mernild and 
others, 2011). 

AIS results 

The mascon solution for the AIS regions is summarized in 
Table 5, and Figures 12 and 13. Unlike the GIS, both the 
magnitude and uncertainties of the GIA corrections for the 
AIS are significant with respect to the total trend. The overall 
AIS trend for the December 2003-December 2010 time 
period studied is — 81±26Gta~ 1 . Our trend estimate 
agrees, within stated uncertainties, with — 69±18Gta _1 
obtained by King and others (2012) for August 2002- 
December 2010, using a comparable recent regional GIA 
model. The AIS trend is dominated by significant mass loss 
from the WAIS, with a trend of —106 ± 1 6 Gta -1 . A trend of 


— 38± 14Gta _1 is found for the AIS peninsula, while the 
EAIS trend is 63 ± 28 Gta -1 . The largest trends are found in 
the WAIS along the Amundsen and Bellingshausen Sea 
coasts (drainage systems 19-22), concentrated at the Pine 
Island embayment (convergence of drainage systems 19- 
22), and at the northern tip of the peninsula (drainage system 
24). These results are in general agreement with —132 ± 60 
Gta -1 for the WAIS and — 60±46Gta~ 1 for the AIS 
peninsula, determined from radar interferometry and re- 
gional climate models for the year 2006 (Rignot and others, 
2008). In addition, our AIS peninsula trend is in good 
agreement with — 42±9Gta~’ for January 2003-March 
2009, determined from a GRACE spherical cap mascon 
solution (Ivins and others, 2011). 

The most striking signal, with the greatest potential 
consequence for sea-level rise, is the large accelerated mass 
loss exhibited by the WAIS at —46 ± 6 Gta~ 2 . While we do 
observe interannual variability of the WAIS annual balances, 
there is a consistent negative trend in these data, repre- 
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Fig. 9. The annual mass balance mean (2004-10) and annual anomalies about the mean for GIS (blue), GOA (green), EAIS (purple), WAIS 
(orange) and AIS peninsula (cyan) from the EEMD analysis of the v12 mascon solution. The values and corresponding errors bars are 
determined from the mean extrema and extrema error ellipses shown in Figure 8. 
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Fig. 10. The GIS annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 mascon solution. The 
timing of the seasons are determined by the mean EEMD extrema in Figure 8a, and the mass values are determined by the result of the EEMD 
method (summing the IMFs from the one preceding the 'annual' IMF to the final IMF) for each individual mascon time series at the time of 
the total GIS extrema. 


senting a persistent multi-year acceleration of mass loss 
(Figs 8c! and 9). The largest accelerated loss is concentrated 
at Pine Island embayment and along the Amundsen and 
Bellingshausen Sea coasts with significant mass loss 
extending to the Zumberge coast (drainage system 1) during 
the 2009 balance year (Fig. 14). These changes are consist- 
ent with observations of grounding-line retreat (Rignot and 
others, 2008), acceleration (Joughin and others, 2010) and 
drawdown (Wingham and others, 2009) of the Pine Island 
Glacier driven by basal melting of its large ice shelf 
(Pritchard and others, 2012). 

The EAIS trend (63±28Gta _1 ) and accelerated mass 
gain (22 ± 22 Gta~ 2 ) (Table 5 and Fig. 8c), are not persistent 
multi-year signals, but are due to an exceptional positive 
surface mass-balance anomaly during the 2009 accumu- 
lation season, concentrated primarily along the coast of 
Queen Maud Land (drainage systems 5-7), but also along 
the coast of Wilkes Land (drainage systems 11 and 12) 
(Figs 8c and 14). These observations have similar patterns to 
melting days anomalies derived from passive microwave 
observations which, for the entire AIS, reached a new 
historical minimum for the 1980-2009 period (Tedesco and 
Monaghan, 2009). In addition, our GRACE observations 
agree quite well with the timing and spatial location of a 
significant 2009 accumulation anomaly observed in the 
surface mass-balance record derived from the regional 
atmospheric climate model RACM02 (Lenaerts and others, 


2012; personal communication from J.T.M. Lenaerts, 
[AUTHOR: Please provide date]). 

The northern tip of the AIS peninsula exhibits consistent 
negative annual balance, with the exception of the 2008 
balance year (Fig. 14), which gives rise to the shallow 
positive acceleration shown in Figure 13. The southern area 
of the peninsula shows small positive and negative annual 
balances, with the exception of large negative annual 
balances for 2007 and 2009, giving rise to the sharp losses 
shown in the AIS peninsula time series for these years 
(Fig. 8e). These large negative annual balance anomalies 
give rise to the acceleration signal in the southern peninsula 
shown in Figure 13, and the overall peninsula acceleration 
of — 9±5Gta~ 2 . Large mass losses on the AIS peninsula 
result from tributary glaciers that feed into ice shelves that 
have recently disintegrated (Shuman and others, 201 1), and 
increasingly negative mass balances resulting from strong 
climatic warming in the region. Our observations agree with 
satellite observations of surface lowering (Fricker and 
Padman, 2012) and model estimates of increasing mass 
losses from peripheral glaciers and ice caps in the region 
(Flock and others, 2009). 

GOA results 

The mascon solution for the GOA glacier region is 
summarized in Table 5 and Figures 8b and 15. The overall 
trend of the GOA glaciers for the time period studied is 
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Fig. 11. The v12 10 day mascon solution time series (circles) and 
10 day Gaussian smoothed solution (curves) for selected drainage 
systems within the GIS and AIS, and mascons within the GOA land- 
ice region. The colors of the plotted curves correspond to the 
drainage system and mascon colors in Figure 1. 


—69 ±11 Gta -1 , and —77 ± 1 1 Gta -1 if we take the mean 
of the 2004-10 balance years (Fig. 9). Most of the GIA 
correction in this region results from uplift corrections clue to 
LI A mass loss in the Glacier Bay region of Alaska. Our 
estimates very nearly agree to within stated error bars with 
— 46±7Gta~ 1 reported by Jacob and others (2012). 
Differences in the best estimates between these two studies 
likely results from different treatment of TWS, and slightly 
different solution domains. 

GOA glaciers exhibited large seasonal and interannual 
variability during the 2003-1 0 measurement period (Figs 8b, 
16 and 11). We observed large mass losses during summer 
2004, in response to record warm temperatures across 
Alaska that year (Truffer and others, 2005). Even larger mass 
losses occurred in 2009, despite the absence of a strong 
temperature forcing, and we attribute these losses to the 31 
March 2009 eruption of Mt Redoubt that spread volcanic 
ash across much of the GOA region (Schaefer and others, 
2012) and likely enhanced melt rates through a reduction in 
surface albedo (Arendt and others, submittedfAUTHOR: 
Only accepted papers can be cited. See note in reference 
list]). Large 2007-08 winter accumulation was followed by a 
cool 2008 summer season, resulting in near-balance 
conditions during 2008. This offset the strong 2004 and 
2009 melt seasons and resulted in no significant accelera- 
tion in GOA mass balance during the GRACE period. 

The spatial distribution of mass balances largely follows 
the distribution of glacierization in the GOA region, with 
largest mass losses occurring in the St Elias and Glacier Bay 
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Fig. 12. The v12 10 day mascon solution time series (circles) and 
10 day Gaussian smoothed solution (curves) for the AIS total (blue), 
AIS <2000 m elevation (purple), AIS >2000 m elevation (green), 
EAIS (orange), WAIS (cyan) and AIS peninsula (gold). The red curves 
show the best-fit trend and acceleration. 
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Fig. 13. The AIS v12 mascon solution mass changes computed over 
the time period 1 December 2003-1 December 2010. (a) Trend 
corrected with the IJ05 R2 GIA model and (b) acceleration. The 
EAIS, WAIS and AIS peninsula regions are demarcated by the cyan 
borders. 
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Fig. 14. The AIS annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 mascon solution. The 
timing of the seasons are determined by the mean EEMD extrema in Figure 8c, d and e, and the mass values are determined by the result of 
the EEMD method (summing the IMFs from the one preceding the 'annual' IMF to the final IMF) for each individual mascon time series at the 
time of the extrema for the region in which the mascon is located (east, west or peninsula). The EAIS, WAIS and AIS peninsula regions are 
demarcated by the cyan borders. The mean annual mass balance is shown for the v12 solution corrected with the IJ05 R2 GIA model. 
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regions. The highly variable distribution of glaciers in the 
GOA region, with large signals occurring in the center of the 
domain and much smaller signals at the northwestern and 
southeastern edges, likely results in smearing and leakage of 
signal between individual mascons. Therefore, although the 
total GOA signal is not affected by these problems, 
additional processing steps are required to further constrain 
mass changes over drainage systems smaller than our 
current solution domain (Arendt and others, submitte- 
d[AUTHOR: Only accepted papers can be cited. See note 
in reference list]). 


Fig. 15. The GOA mascon trends computed from the v12 mascon 
solution corrected with the ICE5G_Paulson GIA model and the 
Larsen LIA model over the time period 1 December 2003- 
1 December 2010. 


CONCLUDING REMARKS 

We have developed a new global mascon solution spe- 
cifically tuned to observe high-latitude land-ice mass 
evolution estimated directly from the GRACE intersatellite 




cm w.e. 


Fig. 16. The GOA land-ice region annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 
mascon solution. The timing of the seasons is determined by the mean EEMD extrema in Figure 8b, and the mass values are determined by 
the result of the EEMD method (summing the IMFs from the one preceding the 'annual' IMF to the final IMF) for each individual mascon time 
series at the time of the total GOA extrema. 


3B2 v8.07jA/V 12th June 2013 Article ref 1 2J 147 Typeset by Sukie Proof no 1 


Luthcke and others: Land-ice evolution from GRACE 


17 


range-rate data. The study is unique in many ways including: 
(1) a detailed error analysis of a constrained global mascon 
solution with estimates of signal leakage in and out of the 
regions of interest; (2) a controlled comparison to solutions 
derived from monthly spherical harmonic solutions using 
the same GRACE data processing and forward models, 
providing estimates of uncertainty due to solution technique; 
(3) fully iterating a global mascon solution to convergence 
and quantifying the signal recovery improvement and 
GRACE intersatellite range-rate residual reduction and (4) 
applying the EEMD to compute land-ice annual balances. 
We estimate the land-ice mass trends for the time period 
1 December 2003-1 December 2010 to be: —230 ± 12 Gt 
a -1 for the GIS; —81 ± 26 Gta~’ for the AIS and —69 ± 1 1 
Gta -1 for the GOA glaciers. For the subregions of the AIS we 
estimate: 63±28Gta~ 1 for the EAIS; — 1 06 ± 1 6 Gta -1 for 
the WAIS and — 38±14Gta~ 1 for the AIS peninsula. The 
total over the land-ice regions studied is —380 ±31 Gta~\ 
equivalent to — 1 .05 ± 0.09 mm a -1 sea-level rise. Over the 
same time period we estimate the mass acceleration to be: 
— 1 0 ± 6 Gta~ 2 for the GIS; —33 ± 26 Gta~ 2 for the AIS and 
2±7Gta~ 2 for the GOA glaciers, giving a total of 
—41 ±27Gta~ 2 , equivalent to —0.11 ±0.08 mm a -2 sea- 
level rise. For the subregions of the AIS we estimate the 
accelerations to be: 22±22Gta~ 2 for the EAIS, — 46 ± 6 
Gta~ 2 for the WAIS and — 9 ± 5 Gta~ 2 for the AIS peninsula. 
The trends and accelerations determined over our study 
period have been shown to be highly dependent on 
significant seasonal and annual balance anomalies, making 
it difficult to predict future land-ice mass balance and its 
contribution to sea level. One possible exception is the 
WAIS acceleration signal, which shows a persistent negative 
trend in annual balance. The iterated global mascon solution 
presented here provides the spatial and temporal monitoring 
of land-ice evolution to aid in improved understanding and 
modeling. 
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